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ABSTRACT 

Wc propose the multi-frequency synthesis (MFS) algorithm with spectral correction 
of frequency-dependent source brightness distribution based on maximum entropy 
method. In order to take into account the spectral terms of n-th order in the Taylor 
expansion for the frequency-dependent brightness distribution, we use a generalized 
form of the maximum entropy method suitable for reconstruction of not only positive- 
definite functions, but also sign-variable ones. The proposed algorithm is aimed at 
producing both improved total intensity image and two-dimensional spectral index 
distribution over the source. We consider also the problem of frequency-dependent 
variation of the radio core positions of self-absorbed active galactic nuclei, which should 
be taken into account in a correct multi-frequency synthesis. First, the proposed MFS 
algorithm has been tested on simulated data and then applied to four-frequency syn- 
thesis imaging of the radio source 09544-658 from VLBA observational data obtained 
quasi-simultaneously at 5, 8, 15 and 22 GHz. 

Key words: Methods: data analysis; techniques: interferometric image processing, 
high angular resolution; galaxies: nuclei, jets. 



1 INTRODUCTION 

At present, very-long-baseline interferometry (VLBI) is the 
most powerful tool for studying the morphological structures 
as well as kinematic, polarization, and spectral characteris- 
tics of active galactic nuclei (AGN). It allows objects to be 
imaged with a very high angular resolution reaching frac- 
tions of a milliarcsecond (mas). One of the topical problems 
of VLBI mapping of AGN is multi-frequency image synthe- 
sis. Our interest in this method is mainly related to the 
peculiar geometry of the future high-orbit ground-space ra- 
dio interferometer Radioastron (Kardaslicv I QQtI ), which is 
expected to provide an ultra-hi gh resolution (microarcsec- 
onds), but poor aperture filling (|Baikovall2005l ). 

Multi-frequency synthesis in VLBI suggests mapping 
AGN at several frequencies simultaneously to improve the 
instrument aperture filling. This is possible, because the in- 
terferometer baselines are measured in wavelengths of the 
emission being received. The problem of multi-frequency 
synthesis is complicated due to the frequency-dependent 
brightness of a radio source. Hence, to avoid undesirable ar- 
tifacts in the reconstructed image, spectral correction should 
be made at the stage of its deconvolution. 
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ISault fc Wieringal (| 19941 ). and ISault fc Conwavl (|l999l ) m- 
vestigated the influence of spectral effects on the image and 
developed the methods of their correction. These authors 
showed that if narrow frequency bands, up to ±12.5% of 
the reference frequency, are used, the effects of the spectral 
dependence of the brightness of a radio source can usually 
be ignored for dynamic ranges of less than 1000:1. When the 
spectral errors are above the noise they can be recognized 
and removed to ensure required dynamical range of images. 
The spectral errors can usually be accounted for by param- 
eterizing the image in terms of two unknowns at each pixel: 
the intensity at some reference frequency and spectral index 
if spectral variation of the source emission is modelled by a 
power-law relationship. Thus, the use of MFS doubles the 
number of unknowns but in case of MFS we have n/ more 
data of observations, where rtf i s num ber of frequencies. As 
discussed in paper bv I Conwavl (|l99ll ). if n/ greater than 2 
the MFS image remains better constrained than the single 
frequency image. 

The algorith m of linear spec tral correction based on the 
CLEAN method jHogb omlll974) and called "do uble decon- 
volution" (|Conwav. Cornwell fcWilkinsonll 19901 ) is the best- 
studied one. In this algorithm, the "dirty" image is first de- 
convolved with an ordinary "dirty" beam and the residual 
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map is then deconvolved with the beam responsible for the 
first-order spectral term. An improvement of this method 
consisting in simultaneous reconstruction of the sought-for 
image and the map of t he spectral term was proposed by 
ISault fc Wieringal (Il994 ). The vector relaxation a lgorithm 
developed bv iLikhachev. Ladvnin fc GuirinI (|2006l ) may be 
considered as a generalized CLEAN deconvolution method 
that can account for spectral terms of any order. Description 
of the MP S technique i n a co mmo n ma t hema tical framework 
is given in lRau et all (|2009l ) and lRaul (|201Ct l. 

An alternative deconvolution method also actively 
used in radio astronomy is the maxim um entr o py m ethod 
(MEM). ME M was first proposed by iFriedenI j 19721) and 
lAblesI (| 19741 ') for the reconstruction of images in optics and 
radio astronomy, respectively. Sinc e then the method has 
been developed in many studies (ISkiUing fc Brvan' 'l984l: 
Naravan fc Nitvananda 1986.: Wcrncckc fc D'Addario 1976 : 



CornweU fc Evanslll985l : iFrieden fc Baikoval 119941 ) and im 
plemented in a number of software packages designed for 
image reconstruction (MEMSYS, AIPS, etc.). A compar- 
ative analysis of CLEAN and M EM in VLBI is given in 
ICornwell. Braun fc Brigg^ (|l999l ) ; essentially, these methods 
complement each other. In particular, CLEAN is preferred 
for reconstructing images of compact sources from relatively 
poor data, while MEM is more suitable for imaging of ex- 
tended sources from better-quality data. 

A severe drawba ck of MEM compared to CLEAN , the 
bias of the solution (jCornwell. Braun fc Brigg^ Il999l ) , can 
easily be removed by generalizing the met hod to enable 
the reconstruction of sign -variable functions (|Baikovall 19921 : 
iFrieden fc Baikoval[l99i ). This generalized MEM also per- 
mits difference imaging making it possible to substantially 
broaden the dynamic range of maps of sources, i ncluding 
both compact and extended , faint components (jBaikoval 
I2OO7I : iRastorgueva et aLll201ll). 

As shown by Baikoval ( 20081 ) . applying Shannon's max- 
imum entropy method allows a tangible progress to be 
achieved in solving the problem of multi-frequency synthesis 
owing to the possibility of a simple allowance for the spectral 
terms of any order. This, in turn, allows the range of synthe- 
sized frequencies to be extended significantly. However, the 
multi-frequency synthesis algorithms based on both CLEAN 
and MEM deconvolution discussed here can be directly ap- 
plied only to those radio sources for which no frequency- 
dependent image shift is observed. Otherwise, as shown by 
ICroke fc Gabuzdal (|2008l ). an additional operation to align 
images at different frequencies should be performed to ob- 
tain the proper results in a multi-frequency data analysis. 

In this paper, our goal is to deduce our multi-frequency 
synthesis algorithm based on MEM and to show the im- 
portance of the procedure for precorrecting the frequency- 
dependent image shift while implementing multi-frequency 
synthesis. 

The paper is structured as follows. The frequency de- 
pendence of the image of a radio source is described in 
Section (2] Frequency-dependent constraints on the visibil- 
ity function are derived in Section [3] The reconstruction 
method is given in Section O Prior to describing the MFS 
algorithm, we consider the maximum entropy method and 
its generalized form in Subsections 14.11 and 14.21 respectively. 
The proposed multi-frequency synthesis algorithm with fre- 
quency correction is deduced in Subsection 14.31 We test our 



algorithm in a number of model experiments in Section (5) 
Discussion of the problem of aligning frequency-dependent 
images to properly construct the spectral index distribution 
is given in Section |6l And, finally, in Section [T] we present 
the results of applying the MFS algorithm proposed to four- 
frequency VLBA data for BL Lacertae object 0954-1-658. 



2 FREQUENCY DEPENDENCE OF THE AGN 
RADIO BRIGHTNESS 

The dependence of the intensity of a radio source on fre- 
quency v in the model of synchrotron radiation is given by 
(|Conwav. CornweU fc Wilkinsonlll990l ) 



liu) = liuo) - 



(1) 



where /(i^o) is the intensity of the radiation at the reference 
frequency uo and a is the spectral index. To simplify the 
writing, we will hereafter set 7o — liyo). 

Retaining the first Q terms in the Taylor expansion 
of ([1} at point vq, we can write the following approximate 
equality: 



9 = 1 



(2) 



where 



a(a-l)---[a-(g-l)] 



From ([2]), for each pixel (k,l) of the source's two- 
dimensional (A^ X A'^) brightness distribution we have 



i{k,i)^io{k,i) + Y,i,{k,i) 



9=1 



(3) 



where k,l = 1, . . . , N . 

Thus, the derived brightness distribution over the 
source Q is the sum of the brightness distribution at the 
reference frequency vq and the spectral terms with the gth- 
order spectral map depending on the spectral index distri- 
bution over the source as follows: 

J,(fc,0=Jo(fc,0 "(^'^^---["(^,'^^'(^'^^' . (4) 

Of greatest interest is the first-order spectral map 

Ii{k,l) = Io{k,l)a{k,l), 

because the spectral index distribution over the source can 
be estimated from that in the following way: 

a{k,l)^h{k,l)/h{k,l). (5) 



3 CONSTRAINTS ON THE VISIBILITY 
FUNCTION 

The complex visibility function is the Fourier transform of 
the intensity distribution over the source that satisfies the 
spectral dependence ^ at each pixel of the map (fc, I). Given 
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the finite number of terms in the Taylor expansion ((3]), the 
constraints on the visibiHty function V can be written as 

=F{/(fc,/)} X D„,,,„ (6) 

where F denotes the Fourier transform and D represents the 
transfer function, which is the 5-function of it and v for each 
measurement of the visibility function; different sets of 5- 
functions correspond to different frequencies v, as suggested 
by the indices of u and v. 

Let us rewrite (|6]) for the real and imaginary parts of 
the visibility function Ki„,ti„ = + jBu^,v^ by taking 

into account the measurement errors as 

9=0 fc,i ^ ^" ' 

(^^^)'+'?^:,.. = B..,.., (8) 

q=0 k,l \ '^0 J 

where and 6u™,„„ are the constant coefficients (cosines 

and sines) that correspond to the Fourier transform, rf^^^^^ 
and are the real and imaginary parts of the instru- 

mental additive noise distributed normally with a zero mean 
and a known dispersion Gu^^v^,- 



4 METHOD 

4.1 Mciximum Entropy Method 

MEM is one of a large class of non-linear informational 
methods based on the optimization of a functional speci- 
fied by some informational criterion for the quality of the 
solution subject to the constraints that flow from the data. 
In our case, maximizing the Shannon entropy consists in 
finding the maximum of the functional 

E = - j x{t)\n[x{t)) At, (9) 

where x{t) is the desired distribution. 

Since imaging in VLBI implies dealing with digital 
data, we present a discrete formulation of the optimiza- 
tion. Let a map of an object with a finite carrier be 
sampled in accordance with Kotelnikov-Nyquist theorem 
jOppenheim fc Schafer|[l999l ) and have a size oi N x N pix- 
els. We denote the discrete measurements of the desired dis- 
tribution by 

Xki, k,l = 1, . . . ,N - 1. 

We denote the known measurements of the two-dimensional 
Fourier spectrum of the object, which represent the visibil- 
ity data, in accordance with the Van Cittert-Cernike theo- 
rem, as follows, separating the real. Am, and imaginary, Bm, 
parts: 

Vm^Am+jBm, m=l,...,M, 

where M is the number of known measurements and m is 
the number of the current measurement with coordinates 
(umjVm.) in the uv-plane, not necessarily located at nodes of 



the coordinate grid. This last circumstance means that there 
is no problem with pixelization of the data in the frequency 
domain, which represents a certain technical advantage of 
this method over other methods and appreciably enhances 
the accuracy of the reconstruction. 

The practical MEM algorit hm we applied , taking into 
account the errors in the data (|Baikovalll993l ). implies the 
solution of the conditional optimization problem 

min^^x.,h.(x-.0 +^_Y^ (v-.r + (r>l^)\ (10) 

k I m 

'^'^Xkia'^i +ri" = A,n, (11) 

k I 

^^SfcifellJ +7?™ = (12) 

k I 

Xki > 0. (13) 

As we can see from (jlOp . the optimized functional has 
two parts: a Shannon entropy functional and a functional 
that is an estimate of the difference between the recon- 
structed spectrum and the measured data according to the 
criterion. This latter functional can be considered an ad- 
ditional regulating, or stabilizing, term acting to provide a 
further regularization of the MEM solution. The influence of 
this additional term on the resolution of the reconstruction 
algorithm must be kept in mind. 

Equations pil) . (|12p represent linear constraints on the 
unknown images Xki as well as noise terms rj^ and r;™. The 
non-negativity constraint p3[) on the image can be omitted 
in this case due to the nature of the entropy solution, which 
is purely positive. If the total flux of the source Fo is known, 
this automatically leads to the normalization of the solution: 

k I 

The numerical algorithm for solution pop - (|12|> . treated as a 
non-linear optimization problem based on t he method of La - 
grange multipliers, is considered in detail in lBaikoval(|200'i1) . 
Here, we present only the solution: 

Xki = exp{-'^{amaTi + Pnib]^i) - 1), (14) 

m 

rc 2 im ^ R 

expressed in terms of the Lagrange multipliers (dual vari- 
ables) Um and /3„i, through which the constraints (jlip and 
(|12|l . respectively, enter the Lagrange functional. 

As we can see from (|14|l . the standard MEM image is 
evidently positive. It can be shown that the MEM Hesse ma- 
trices everywhere are positive-definite, so that the entropy 
functional is convex and the solution is global. Various gra- 
dient methods can be used to search for the extrema of the 
corresponding dual functional. We use a coordinate-descent 
method. 

4.2 Generalized Msixinium Entropy Method 

The GMEM was designed for the re construct i on of 
sign-variable and complex functions (jBaikoval Il992l : 
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iFrieden fc Baikoval [l993 : iBaikoval l2007l ). For the GMEM, 
dealing with sign- variable real distributions, the Shannon- 
entropy functional has the form 

E = - j {x'^{t) \n{[ax:^ it)] + {t)\n[ax^ it)]} At, 

where x'^{t) and x^ [t] ^ are the positive and negative 
components of the sought-for image x{t), i.e. the equation 
x{t) = x'^ (t) — x^ (t) holds; a > is a parameter responsible 
for the accuracy of the separation of the negative and pos- 
itive components of the solution x{t), and therefor e critical 
for th e resulting image fidelity. As it was shown bv lBaikoval 
(|2007l ) solutions for x'^{t) and x (t) obtained with the La- 
grange optimization method are connected by the expression 

x+{t) ■ x~{t) = exp(-2 - 21na) = K{a), 

which depends only on the parameter a. 

This parameter is responsible for dividing the positive 
and negative parts of the solution: the larger a allows the 
more accurate discrimination (since K{a) — as a — s- oo). 
On the other hand, the value of a is constrained by com- 
putational limitations. The main constraint comes from the 
term in the optimized functional, which depends on the 
data errors. The larger a standard deviation is, the higher 
value of a could be set. If data are very accurate, a lower 
value of a is needed. In practice, a is chosen empirically. In 
our case we had to compromise between data errors, which 
determine resolution of the final MEM solution, and a need 
to divide the positive and negative parts of the solution as 
accurately as possible. It is fair to say that given fixed errors 
in the data, a maximum possible chosen value of a provides 
us with the best possible resolution of the MEM-solution. 
In this work we used a = 10®. 



4.3 GMEM-based MPS algorithm 

In this case, the distributions Iq{k,l), q — 0, . . . ,Q — 1; 
k,l — 1, . . . , N, and the measurement errors of the visibil- 
ity function r]'j^,^ v^^, are unknown. Note that although 
the brightness distribution over the source is described by a 
non-negative function, the spectral maps of arbitrary order 
((U can generally take both positive and negative values be- 
cause the spectral index distribution over the source is an 
alternating one. 

By setting, in accordance with the approach described 
above: 

I) = I+{k, l)-I-{k,l), g = 1, . . . , Q - 1, 
we obtain the following entropic functional to be minimized: 

E = ^7o(A:,01n[a/o(fc,0] (15) 



Table 1 . Parameters of model source images at different frequen- 
cies. 



+ Yi^ik, I) ln[a/+(fc, 0] + I^{k, I) ln[a/-(fc, /)]} 

9=1 k,l 



/o(fc,0 ^ 0, /+(fc,0 > 0, I~ik,l) > 0. 



Frequency 


5'tot 


'S'peak 


Entropy 


(GHz) 


(Jy) 


(Jy/pixel) 




5.0 


5.17 


0.187 


-17.8 


8.4 


4.76 


0.192 


-16.2 


13.6 


4.64 


0.196 


-15.5 


15.3 


4.63 


0.197 


-15.4 


22.2 


4.65 


0.201 


-15.3 



The linear constraints ([7]) and ([8]) on the measured vis- 
ibility function will be rewritten accordingly: 



y^Jo{k,l)atl,v^ 

k.l 



(16) 



Q-i 



9=1 k.l 

V — J/o \ ' re _ - 

fO J 

^ Jo(fc,o&L",„„ 

k.l 

Q-1 



(17) 



9=1 k.l 

fl^-niY , im o 
\ J 

Minimizing the functional (|15[) with constraints p6p - 
(|17[) constitutes the essence of the MEM-based multi- 
frequency synthesis algorithm that seeks the solution for 
all unknown Jo (fc, 0) lt'^~\k,l), 9 = 1)---,Q — 1, k,l = 
1, . . . ,N and ffu^^v^ , '7u™,d„ • A. detailed algorithm for numer- 
ical implementation of the proposed m ulti-frequency synthe- 
sis method is given in iBaikoval (|2008l ). 



TESTING THE METHOD. 
RESULTS 



SIMULATION 



Here we present the results of testing our MFS deconvolution 
algorithm on the example of four-frequency synthesis using 
simulated VLBI data at four frequencies of 5, 8, 15 and 
22 GHz. As a reference frequency in the MFS algorithm, we 
adopted the central frequency equal to 13.6 GHz. A model 
source map at this frequency and a model spectral index 
distribution over the source are presented in Fig. [TJ left and 
right respectively. 

As one can see, the source shows a structure on a mil- 
liarcsecond angular scale consisting of a bright compact core 
and a one-sided jet. Note that the model map contains also 
a number of weak small-scale details scattered around the 
main structure. Such a complication of the source structure 
was made in order to test the ultimate capabilities of the 
MFS algorithm. Note also that the structure of the model 
source was built similar to the structure of the radio source 
0954-f658 which will be considered in Section [T] 
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The Mj>plane coverages related to a model interferom- 
eter consisting of 10 baselines of the VLBA array (BR-LA, 
BR-MK, BR-NL, BR-PT, BR-OV, BR-SC, BR-KP, BR-HN, 
BR-FD, LA-MK) and source with a declination S = 70° for 
each of "observation" frequencies are shown in Fig. [2] Note 
that the choice of baselines was not principal for our task 
and was made in an arbitrary way. It was assumed that a 
duration of the observations is nine hours and visibility data 
are formed every 30 minutes. As it can be seen, the uv-p\a.ne 
at different frequencies has the same topology but scaled 
accordingly. As far as interferometer baselines are measured 
in wavelengths, u- and v- visibility coordinates are propor- 
tional to an observation frequency. 

Model maps of the source at frequencies of 5, 8, 15 and 
22 GHz are shown in Fig. [3l left column. These maps were 
obtained from the model map shown in Fig. [1] left, accord- 
ing to the model spectral index distribution (Fig. [T] right) 
and expression ^ for the spectral dependence. Parameters 
of the model source intensity maps, such as total flux density 
(Stot), peak flux density (Speak) and entropy, are given in Ta- 
ble [1] Comparison of the maps and their parameters shows 
although not strong but still quite noticeable frequency de- 
pendence. 

The complex visibility functions, computed in accor- 
dance with the Van Cittert-Cernike theorem, consists of 
172 samples for each "observation" frequency. Each visibil- 
ity value was aggravated with additive random error to form 
data with a typical signal-to-noise ratio of averaged and self- 
calibrated VLBI data of about lO**. It is necessary to empha- 
size that we do not consider here the se lfcalibration problems 
IConwav. Cornwell fc Wilkinsonlll990l ) and concentrate only 
on removing spectral errors. 

The main goals of the simulation fulfilled were: (i) to 
show efficiency of the multi-frequency synthesis for improv- 
ing intensity images in case of small-element interferometers 
with poor iii»-plane filling; (ii) to illustrate the consequences 
of multi- frequency synthesis without any spectral correction; 
(iii) to estimate the possibility of reconstruction of spectral 
index maps with a satisfactory quality. 

We performed three tests. The first one was single- 
frequency synthesis of the source images at 5, 8, 15 and 
22 GHz. The second one was multi(four)-frequency synthe- 
sis without any spectral correction. And, finally, the third 
experiment was devoted directly to multi-frequency synthe- 
sis with spectral correction, at the reference frequency of 
13.6 GHz. 

The single- frequency maps are presented in Fig. O right 
column. Their parameters are given in Tabled Analysis of 
the results obtained shows the following. The amount of data 
proves to be too small, and ^i^^plane filling too poor to obtain 
images with sufficient quality. Comparison of Table [T] and [2] 
allows us to judge about distortions of the reconstructed 
images. As a criterium of the reconstructed image quality, 
we chose the signal-to-noise ratio (SNR) listed in the last 
column in Table [2j which was calculated in the following 
way: 




SNR = 



Figure 1. Intensity map for the model radio source at 13.6 GHz 
(left) and model spectral index map (right). Relative right ascen- 
sion and declination are given in mas along the horizontal and 
vertical axes, respectively. Contour levels: 0.2, 0.4, 0.8, 1.6, 3.2, 
6.4, 12.8, 25.6, 51.2, 90% of the peak flux density of 0.196 Jy/pixel. 



where /mod is the model map, /rec is the reconstructed map, 
k,l = 1, . . . , N , where N is the linear size of the map. 

As it is seen from the single-frequency maps, at lower 
observation frequencies, the lower resolution is ensured. The 
lower-frequency maps show a larger-scale structure, while 
the higher-frequency maps show smaller-scale features. The 
highest accuracy of reconstruction is achieved at frequencies 
of 8 and 15 GHz (SNR ~ 16). At the lowest frequency, 
5 GHz, and the highest frequency, 22 GHz, the obtained 
reconstruction quality was much worse (SNR is about 8). 

The iii»-plane corresponding to multi(four)-frequency 
synthesis is presented in Fig. |4l left. The multi-frequency 
map obtained without any spectral corrections of the 
frequency-dependent source brightness distribution is shown 
in Fig.|4l right. Parameters of the synthesized map are given 
in Table [3] Comparing with the model map (Fig.[T}, we can 
see large image distortions which are reflected in such map 
parameters as the total flux density, entropy and SNR (~ 
7). Large image distortions in the case of relatively small 
values of spectral index distribution (Fig. [T] right) can be 
explained by a wide frequency range of the data. 

The intensity map and spectral index image obtained 
using the multi-frequency synthesis algorithm with spectral 
correction are shown in Fig. [5] Note that we tried the dif- 
ferent numbers of spectral terms in expansion ([2]). Here we 
present the results related to utilization of three spectral 
terms. Further increase in the number of spectral terms 
did not improve our results substantially. Using fewer spec- 
tral terms proved to be insufficient due to the wide fre- 
quency range of the data. Having analyzed the results (Fig. [5] 
and Table |3}, we can conclude that taking into account 
the frequency dependence of the source brightness distri- 
bution allowed us to obtain both the source intensity map 
(SNR ~ 36) and spectral index distribution over the source 
with high accuracy. 

Thus, in the source model with sufficiently complicated 
extended structure, typical for AGN at milliarcsecond scales, 
with typical spectral index values, we demonstrated the abil- 
ity of the MEM-based multi-frequency synthesis algorithm 
with spectral correction for both (i) improving intensity im- 
ages and (ii) obtaining spectral index maps of high quality. 
We emphasize that the use of the MFS algorithm is espe- 
cially effective in the case of small-element interferometers 
with poor uv-plane filling. We also showed the consequences 
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Figure 2. Coverages of the uv-plane of the simulated 10-baseline 
interferometer related to data formed at 5 (a), 8 (b), 15 (c) and 
22 (d) GHz. The axes are in units of miUions of wavelengths. 



Table 2. Parameters of single- frequency reconstructed images. 



Frequency 


Stot 


'5'pcak 


Entropy 


SNR 


(GHz) 


(Jy) 


(Jy/pixel) 






5.0 


5.16 


0.168 


-18.3 


7.4 


8.4 


4.71 


0.187 


-16.1 


16.8 


15.4 


4.54 


0.197 


-15.5 


16.5 


22.2 


4.41 


0.189 


-15.7 


8.4 



of ignoring frequency dependence of the source brightness 
distribution in mult i- frequency synthesis algorithm. 



6 THE PROBLEM OF IMAGE ALIGNMENT 

One of the most important sources of information about 
the physical conditions in the radio-emitting regions of 
AGN is the spectral index distribution over the source. 
The core region is usually characterized by a large opti- 
cal depth and an almost flat or inverted spectrum, while 
the jets are optically thin with re spect to synchrotron ra- 
diation and have steeper spectra jCroke fc Gabuz dal l2008l : 
lO'SuUivan fc Gabuzdall2009l : IPushkarev et al.ll2005l ). 



The spectral index distribution over the source can be 
constructed by various methods. The traditional method 

Table 3. Parameters of four-frequency reconstructed images. 



MFS experiment 


•Stot 


'5'pcak 


Entropy 


SNR 




(Jy) 


(Jy/pixel) 






No spectral corr. 


5.85 


0.194 


-22.7 


7.2 


With spectral corr. 


4.71 


0.198 


-16.6 


36.2 
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Figure 3. Model (loft column) and reconstructed (right column) 
images at 5 (a), 8 (b), 15 (c) and 22 GHz (d). Relative right 
ascension and declination are given in mas along the horizontal 
and vertical axes, respectively. Contour levels: 0.2, 0.4, 0.8, 1.6, 
3.2, 6.4, 12.8, 25.6, 51.2, 90% of the peak flux density (see Tables 1 
and 2). 



suggests: (i) formation of images at two separate frequencies, 
ui and with the solutions of the deconvolution problem 
(CLEAN or MEM) being convolved with the same clean 
beam corresponding to the lower observation frequency; (ii) 
calculation of the two-dimensional spectral index distribu- 
tion over the source from equation ([1]). Obviously, this se- 
quence of operations is legitimate only when the positions of 
the VLBI cores of sources (not to be confused with the phys- 
ical core of the source that is undetectable due to absorption 
effects) are frequency-independent. 

The image reconstructi on using the iterative selfcali- 
bration procedure is known (|Thompson. Moran fc SwensorJ 
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Figure 4. Coverage of the uii-plane related to multi(four)- 
frequency data (left) and the image reconstructed from the four- 
frequency data ignoring any frequency dependence of the source 
brightness distribution (right). The u- and v-axes are in units of 
millions of wavelengths. Relative right ascension and declination 
are given in mas along the horizontal and vertical axes, respec- 
tively. Contour levels: 0.2, 0.4, 0.8, 1.6, 3.2, 6.4, 12.8, 25.6, 51.2, 
90% of the peak flux density of 0.194 Jy/pixel. 




Figure 5. Intensity image (left) and spectral index distribu- 
tion (right) reconstructed from the four-frequency data using the 
GMEM-based algorithm with spectral correction (reference fre- 
quency is 13.6 GHz). Relative right ascension and declination are 
given in mas along the horizontal and vertical axes, respectively. 
Contour levels: 0.2, 0.4, 0.8, 1.6, 3.2, 6.4, 12.8, 25.6, 51.2, 90% of 
the peak flux density of 0.198 Jy/pixel. 



[2O0I) to lead to the loss of information about the absolute 
position of the source on the sky: during the phase self- 
calibration, the centroid of the object is placed at the phase 
center of the map with coordinates (0,0). However, since 
most of the radio-lo ud AGN are charact e rized by a dom- 
inant compact core faovalev et al.l l2005l : iLister fc HomanI 
I2OO5I : iLee et aLllioOsI ). the VLBI core of the source coin- 
cides with the peak radio brightness of the source in an 
overwhelming majority of cases. 

Nevert heless, the standard theor y of extragalactic ra- 
dio sources (|Blandford fc Koniglll 19791 ) predicts a frequency- 
dependent VLBI core shift due to opacity effects in the 
source's core region. Synchrotron self-absorption takes place 
in an ultra-compact region near the "central engine" of 
AGN, the mechanism of which is most efficient at low fre- 
quencies. As a result, the apparent origin of the jet mani- 
fests itself farther from the physical core along the jet axis 
at lower frequencies (Fig. [6}. This theoretical prediction 
was confirmed by observations: the frequency-dependent 
shi ft in the c ore po sition was measured for several quasars 
by iLobanovl (| 19981 ). In the literature, this phenomenon is 
also actively debated from the viewpoint of the accuracy 




Figure 6. A scheme illustrating the frequency-dependent posi- 
tion shift of the VLBI core (modified from Kovalcv ct al. (2008a)). 



of astrometric measu rements (|Charlotl I2OO2I : iBobolt j I2OO6I : 
iKovalev et~alll2008al ) . 

It thus follows that the multi-frequency data analysis 
must be preceded by the alignment of images at different 
frequencies. This can be achieved in three ways: (i) per- 
forming VLBI observations of the objects under study to- 
gether with reference sources; (ii) finding the parameters of 
the shift of one image relative to the other by aligning com- 
pact optically thin jet features, which are not subject to 
absorption effects to the same extent as in the so urce's core 
jParagi. Feies fc Frevllioool : iKovalev et al.ll2008ah : (in) find- 
ing the shift paraineters using a cross-correlation analysis 
(ICroke fc Gabuzdall2008l ). Being laborious from the view- 
point of performing observations and their subsequent re- 
duction, the first method gives no significant advantage in 
determining the shift and its accuracy; therefore, the second 
and/or third methods are used more often. 

Recall that the alignment procedure implemented by 
shifting one image relative to the other is equivalent to the 
phase correction of the spectrum (or visibility function) of 
the image being shifted relative to the fixed one. The need 
of precorrecting the data for the source's visibility func- 
tion at different frequencies makes the direct use of the 
multi-frequency synthesis algorithm described above prob- 
lematic, because the frequency dependence of the core shift 
is not known in advance. It can be determined by forming 
the images at each frequency and determining the corre- 
sponding shifts. As it was sho wn by Kovalev et a l. (2008b), 
lO'Sullivan fc Gabuzdal (|2009l ) and lSokolovskv et al.l (|201ll ). 
the frequency dependence of the VLBI core position is well 
fitted by a hyperbolic dependence of the form r oc v^^ . Thus, 
our multi-frequency synthesis procedure can be used after 
allowance for the shifts in the positions of the VLBI cores 
at different frequencies and their coordinates relative to 
the phase center and applying the corresponding frequency- 
dependent phase corrections to the visibility function. 



7 REAL DATA PROCESSING. 

FOUR-FREQUENCY IMAGING FOR 
0954+658 

Here we present the results of applying the developed 
multi-frequency image synthesis algorithm to the real 
VLBI data of the extragalactic radio source 0954-1-658, 
a me mber of the complete sample of BL Lacertae ob- 
jects (|Kiihr fc Schmidtlll990l ). Note that 0954+658 is also 
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Table 4. General parameters of 0954+658. 
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Figure 7. Coverages of the «i>-plane at 5 (a), 8 (b), 15 (c) and 
22 (d) GHz for VLBA observations of 0954+658. The axes are in 
units of millions of wavelengths. 



a member of the IFGL catalog of 7-ray bright sources de- 
tected by the Large Area Telescope onboard the Fermi ob- 
servatory and position ally associated wi th the 7-ray source 
IFGL JIOOO. 1+6539 (|Abdo et al.ll2010l '). This source is of 
our interest because it has a typical parsec-scale morphol- 
ogy that includes an optically thick VLBI core and a one- 
sided optically thin jet, which is expected to manifest itself 
in the spectral index distribution. General information on 
this source is given in Table |4] 

The VLBA observations of 0954+658 were carried out 
in a "snapshot" mode in 1997 April (1997.26) simultane- 
ously at four frequencies: 5, 8, 15 and 22 GHz. The data were 
calibrated in the NRAO AIPS package using standard pro- 
cedures. The images were formed within the framework of 
the Pulkovo "VLBIma ger" softwa re package based on a self- 
calibration algorithm (|Cornwell fc Fonialont 1999) in com- 
bination with a GMEM-based deconvolution procedure. 

The MiJ-plane coverages related to the observation fre- 
quencies of 5, 8, 15 and 22 GHz are shown in Fig. [T] The 
MEM-based single- frequency maps are shown in Fig. [H] Pa- 
rameters of these maps, obtained from the MEM-solutions 
by convolution with "clean" beams that determine the sys- 
tem's resolution at each observation frequency, are given in 
Table El 

The parameters of the frequency-dependent image shift 
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Figure 8. Single-frequency maps of the source 0954+658 recon- 
structed from the data at 5 (a), 8 (b), 15 (c) and 22 (d) GHz. 
Relative right ascension and declination are given in mas along 
the horizontal and vertical axes, respectively. The lowest contour 
levels are given in Table [5] for each image, the values of the suc- 
ceeding levels are doubled. 



Table 5. Parameters of single-frequency maps for 0954+658. 
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found by aligning compact features of the optically thin jet 
are given in Table [6l As expected, the direction of the shift 
coincides with the inner jet orientation. The uv-plaxie related 
to the multi(four)-frequency data is shown in Fig. O left. 

First, we processed the multi- frequency data ignoring 
dependence of the source brightness distribution on observa- 
tion frequency. The image obtained is shown in Fig. [5J right. 
Then we applied our MFS algorithm with spectral correction 
to the multi-frequency data, but we did not correct the data 
in accordance with the frequency-dependent VLBI core po- 



Table 6. Frequency-dependent VLBI core position shifting pa- 
rameters with respect to the core position at 22 GHz. 



Frequency, GHz 


Ar, mas 


Position Angle 


15.4 


0.18 
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Figure 9. Coverage of the uii-plane related to multi(four)- 
frequency data (left) and the image reconstructed from the four- 
frequency data ignoring any frequency dependence of the source 
brightness distribution (right). The u- and v-axes are in units of 
millions of wavelengths. Relative right ascension and declination 
are given in mas along the horizontal and vertical axes, respec- 
tively. Total VLBA flux density is 434 mjy. Beam FWHM is 
0.58x0.43 mas, PA = -21?6. Contour levels: 0.2, 0.4, 0.8, 1.6, 
3.2, 6.4, 12.8, 25.6, 51.2, 90% of the peak flux density of 179 
mjy/beam. 
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Figure 10. Intensity (left) and spectral index (right) maps re- 
constructed from the four-frequency data using the GMEM-based 
algorithm with spectral correction, but without preliminary data 
correction for frequency-dependent core position shift (reference 
frequency is 13.6 GHz). Relative right ascension and declination 
are given in mas along the horizontal and vertical axes, respec- 
tively. Total VLBA flux density is 450 mJy. Beam FWHM is 
0.58x0.43 mas, PA = -21?6. Contour levels: 0.2, 0.4, 0.8, 1.6, 
3.2, 6.4, 12.8, 25.6, 51.2, 90% of the peak flux density of 151 
mjy/bcam. 



sition shift found. Results of this multi- frequency synthesis 
obtained at central reference frequency 13.6 GHz are shown 
in Fig. \W\ 

Finally, we applied our multi-frequency synthesis algo- 
rithm with spectral correction to the observational data that 
were first corrected according to the alignment parameters 
(Table [6|. The image obtained at reference frequency of 
13.6 GHz and the corresponding two-dimensional spectral 
index distribution are shown in Fig. 1111 

The intensity images shown in Fig. [9] [10] and [11] are 
obtained by convolving solutions for Jo with appropriate 
gaussian beam. The spectral index images (Fig. \W\ and llip 
are formed by dividing two solutions for Ii and Jo in ac- 
cordance with equation (O after convolving them with the 
same beam. 

It is necessary to note that the main goal of this section 
is to demonstrate the necessity of taking into account the 
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Figure 11. Intensity (left) and spectral index (right) maps recon- 
structed from the four-frequency data using the GMEM-based al- 
gorithm with spectral correction and preliminary data correction 
for frequency-dependent core position shift (reference frequency 
is 13.6 GHz). Relative right ascension and declination are given 
in mas along the horizontal and vertical axes, respectively. Total 
VLBA flux density is 496 mJy. Beam FWHM is 0.58x0.43 mas, 
PA = -21?6. Contour levels: 0.2, 0.4, 0.8, 1.6, 3.2, 6.4, 12.8, 25.6, 
51.2, 90% of the peak flux density of 189 mJy/beam. 



frequency-dependent image shift in order to correctly map 
the spectral index distribution. 

Let us analyze the results obtained. 

The VLBI structure of the radio source 0954H-658 con- 
sists of an optically thick core (Fig. 1111 right) and an opti- 
cally thin jet initially extending in the North- West direction 
up to ~ 2 mas from the VLBI core and then turning to the 
West. As it is seen from single-frequency maps Fig. [8]), the 
apparent extent of the jet reaches about 12 mas and 2 mas 
at the lowest and the highest frequencies, respectively. Total 
flux density of the source varies from 0.6 up to 0.3 Jy at the 
lowest and the highest frequencies, respectively. One can see 
that the images obtained manifest large-scale structure fea- 
tures at lower observation frequencies and small-scale ones 
at higher frequencies. 

Combining the data obtained at different observation 
frequencies allows us, in general, to reconstruct both large- 
and small-scale structures to a larger extent because of a bet- 
ter filled Mil-plane. That was proved, in particular, in Section 
[5] in model experiments. 

The intensity map (Fig. [5| synthesized directly from the 
linearly combined multi-frequency data without any spec- 
tral corrections shows large source structure distortions and 
many artifacts. Multi-frequency synthesis with spectral cor- 
rection, fulfilled first without taking into account frequency- 
dependent VLBI core position shift, shows a slight improve- 
ment of the source intensity map, but it is still distorted, 
especially in the South- West direction (Fig. I10|l . But the 
main consequence of image misalignment is an incorrect re- 
construction of the spectral index distribution. In our case 
we can see that the spectral index map obtained does not 
correspond to the physical meaning of an optically thick core 
and an optically thin jet: there are negative spectral indices 
in the core region and segments with positive spectral in- 
dices in the jet region. From Fig. 1111 we realize that only 
allowance for the real shift found by aligning features of the 
optically thin jet yielded the proper result. As it is seen, 
we managed to reconstruct a more extended jet structure 
(up to 8 mas along the RA axis) than in the case of using 
only the high-frequency data but with the same high angu- 
lar resolution. The spectral index map adequately reflects 
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Figure 12. Slices of tlie spectral index distributions along the 
jet ridge line beginning from the phase center: (1) and (2) arc 
related to the spectral index maps presented in Fig. llOl and Fig. 1111 
respectively; I is the distance along the ridge line. 



index distribution can be mapped with a high quality as 
well. 
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the physical characteristics of the regions of the optically 
thick compact VLBI core and the optically thin extended 
jet. We see a fairly regular structure with smooth transi- 
tions between segments of different intensities along the en- 
tire source. For comparison, one-dimensional slices of the 
spectral index distributions (see Fig. [TOl left and Fig. [TTl 
left) obtained along the jet ridge line are shown in Fig. 1121 
from which we can graphically evaluate undesirable conse- 
quences of ignoring the frequency-dependent VLBI core po- 
sition shift in the MFS algorithm. 



8 CONCLUSIONS 

We developed and tested an efficient multi-frequency image 
synthesis algorithm with the correction for the frequency de- 
pendence of the radio brightness of a source. The algorithm 
is based on the generalized maximum entropy method; it al- 
lows one to take into account the spectral terms of any order 
and to map both a total intensity image and spectral index 
distribution, which is of great importance in investigating 
the physical characteristics of AGN. 

The advantage of the proposed multi-frequency synthe- 
sis algorithm is that the spectral terms of any order can be 
easily taken into account in the entropic functional being 
minimized. This allows for the spectral correction of images 
to be made both in a wide frequency range and for large 
spectral indices. 

We have shown how important the allowance for the 
frequency-dependent image shift is in applying the multi- 
frequency synthesis algorithm. Our conclusions are based on 
the results of processing the multi-frequency VLBA data for 
the BL Lac object 0954+658 with a fairly complex extended 
jet structure, which also manifests itself in the spectral index 
distribution over the source. 

Analysis of the results obtained shows that multi- 
frequency synthesis is an efficient method for improving 
the mapping quality; low-frequency data allow the extended 
structure of a source to be reconstructed more completely, 
while high-frequency data allow a high spatial resolution to 
be achieved. It should also be emphasized that the spectral 
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